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摘要 


本 文 运 用 函数 对 高 层 建筑 每 段 楼 梯 的 下 散 速率 《单位 时 间 内 人 通 
TO 进行 描述 ， 利 用 泛 函 等 相关 知识 和 计算 机 软件 ,给 出 在 指定 时 间 
内 逃生 人 数 、 逃 生效 率 与 各 楼 层 分配 人 数 间 的 关系 ,同时 讨论 玻 散 方 
式 中 风险 带 来 的 扰动 对 下 散 效率 的 影响 。 通 过 建立 这 样 的 模型 ， 试 图 
对 高 层 建 筑 中 人 员 、 上 店面 等 的 分 布 安排 起 一 定 的 指导 作用 。 


Abstract 

This article tries to describe the evacuation efficiency (people 
flowrate per unit of time) in every flight of stairs in a multi-storeyed 
building, using functional knowledge and computer software, and 
analyze the relationship between evacuation efficiency and personnel 
distribution arrangement. At the same time, this article studies how 
the disturbance brought about by the risk in the evacuation process 
influences the quantity of successfully evacuated people in a limited 
period. Building this model, we expect that it can function as, to a 
certain degree, the guidance of personnel distribution arrangement 


in a multi-storeyed building. 
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第 一 章 前 言 


从 2008 年 四 川 汶川 大 地 震 ，2010 年 舟 曲 泥石流 、 上 海 胶州 路 特大 火灾 ， 
到 今年 雅安 地 震 ， 突 发 性 灾害 事件 频 发 ， 带 来 惨重 的 伤亡 。 各 地 纷纷 加 强 高 密度 
场所 如 医院 、 学 校 等 地 的 应 急 下 散 演练 工作 ， 建 筑 物 紧 急 状 态 下 的 稳固 性 也 作为 
基础 设施 建设 时 的 一 个 指标 更 多 地 为 建筑 者 所 重视 。 

但 我 们 发 现 , 为 了 减轻 突 发 灾害 造成 的 人 员 伤 亡 , 学界 往 往 致力 于 在 以 下 两 
方面 进行 优化 : 加 强 建筑 稳固 性 和 优化 疏散 方案 。 因 此 ， 我 们 试图 提出 一 种 新 的 
思路 来 更 好 地 应 对 潜在 的 危险 , 即 通过 调整 各 楼 层 的 人 员 分 配 来 加 快 疏 散 效率 和 
增强 朴 散 过 程 中 的 稳定 性 。 

在 实际 生活 中 ， 高 层 建筑 往往 兼 有 办 公 、 商 贸 、 和 餐饮 等 多 种 功能 ， 而 各 楼 层 
的 功能 直接 影响 了 各 楼 层 正常 状态 下 的 人 口 密 度 , 因此 我 们 可 以 通过 调整 各 楼 层 
的 功能 来 改变 人 员 比 例 。 以 学 校 为 例 ， 可 以 通过 改变 常用 教室 、 多 媒体 教室 、 空 
教室 等 在 各 楼 层 的 分 布 调整 人 员 分 布 , 而 大 型 商城 则 可 以 通过 调整 餐饮 区 、 购 物 
区 、 仓 库 等 的 位 置 实现 调整 。 这 样 的 调整 对 建筑 的 功能 和 经 济 效益 影响 不 大 ,但 
会 较 大 地 改变 紧急 疏散 的 效率 和 安全 性 。 

角 认 了 该 思路 的 现实 可 执行 性 后 , 我 们 着 手 建立 模型 ， 由 最 简单 的 三 层 建筑 
入 手 《〈 第 一 层 视 为 可 以 不 经 楼 梯 直 接 逃 生 )， 尝 试 给 出 描述 各 楼 梯 口 疏散 效率 的 
方程 。 在 建立 方程 的 过 程 中 , 我 们 需要 应 用 到 人 员 平 均 速 度 与 已 有 人 员 数 量 的 关 
A, 查阅 相关 文献 ， 我 们 发 现 这 个 方程 被 通常 地 表示 为 负 指数 函数 的 形式 。 但 是 
经 过 实地 测算 , 我 们 感到 该 方程 并 不 适用 于 我 们 的 模型 ， 因此 对 其 进行 了 优化 改 
R, 采用 更 合理 的 曲线 描述 人 员 平 均 速度 与 已 有 人 员 数 量 的 关系 ,建立 了 方程 组 。 

运用 数学 工具 分 析 , 我 们 得 到 了 人 员 分 配 比 对 总 体 疏 散 效率 的 影响 ， 并 能 够 
寻找 到 实现 最 佳 疏 散 效 率 的 人 员 分 配 比 。 同 时 ,考虑 到 楼 梯 口 处 大 量 人 员 进 入 对 
原 有 玻 散 秩序 可 能 造成 的 “冲击 ”效果 ， 我 们 将 其 用 相应 位 置 的 扰动 进行 描述 ， 
观察 分 析 其 对 最 终 疏 散 效果 的 影响 ， 从 而 得 到 人 员 分配 比 与 疏散 风险 性 的 关联 。 

KA SRB CRS TE, 我 们 可 以 通过 调整 两 者 的 重要 性 权重 , 来 找 
到 最 佳 的 人 员 分 配 比 ， 同 时 满足 两 方面 的 要 求 ， 从 而 指导 建筑 的 内 部 布局 ， 实 现 
对 突 发 灾害 时 人 员 安全 性 的 进一步 保障 。 


zu 


第 二 章 相关 知识 准备 


2.1 偏 导 数 


一 个 多 变量 的 函数 的 偏 导 数 是 它 关 于 其 中 一 个 变量 的 导数 , 而 保持 其 他 变量 
恒定 (相对 于 全 导数 ， 在 其 中 所 有 变量 都 允许 变化 )。 函 数 f 关 于 变量 x A 


数 写 为 六 或 志 。 偏 导数 符号 是 圆 体 字母 5， 区别 于 全 导数 符号 的 正体 d 


2.2 常 微分 方程 与 偏 微分 方程 


常 微分 方程 (Ordinary Differential Equation， 简 称 ODE) 是 未 知 函 数 只 含 
有 一 个 自 变 量 的 微分 方程 。 

偏 微 分 方程 (Partial Differential Equation, 简称 PDE) 指 含 有 未 知 函 数 及 其 
偏 导 数 的 方程 。 描 述 自 变 量 、 未 知 函 数 及 其 偏 导数 之 间 的 关系 。 符 合 这 个 关系 的 
函数 是 方程 的 解 。 偏 微分 方程 分 为 线性 偏 微分 方程 式 与 非 线 性 偏 微分 方程 式 , 常 
常 有 几 个 解 而 且 涉 及 额外 的 边界 条 件 。 


2.3 变 分 法 


变 分 法 是 处 理 泛 函 的 数学 领域 ,， 和 处 理 函 数 的 普通 微 积分 相对 。 壁 如， 这 样 
的 泛 函 可 以 通过 未 知 函 数 的 积分 和 它 的 导数 来 构造 。 变 分 法 最 终 寻求 的 是 极 值 函 
数 : 它们 使 得 泛 函 取得 极 大 或 极 小 值 。 


2.4 线性 代数 


线性 Cinear) 指 量 与 量 之 间 按 比例 、 成 直线 的 关系 ， 在 数学 上 可 以 理解 为 
一 阶 导 数 为 常数 的 函数 。 非 线性 (non-linear) 则 指 不 按 比 例 、 不 成 直线 的 关系 ， 
一 阶 导数 不 为 常数 。 

线性 代数 起 源 于 对 二 维和 三 维 直 角 坐 标 系 的 研究 。 在 这 里 ， 一 个 向 量 是 一 个 
有 方向 的 线段 ,由 长 度 和 方向 同时 表示 。 ALA AE ACR EI SBR ER 
限 维 空间 。 一 个 维 数 为 n 的 向 量 空 间 叫 做 n 维 空间 。 在 二 维和 三 维 空间 中 大 
多 数 有 用 的 结论 可 以 扩展 到 这 些 高 维 空间 。 


ig 


第 三 章 建立 模型 
3.1 NEES Leet BA 


3.1.1 函数 提出 过 程 


在 通道 不 太 长 以 至 于 涌 入 人 员 对 通道 流量 的 滞后 效应 不 明显 时 , 可 以 视 通 道 
流量 与 人 员 密 度 成 确定 的 函数 关系 。 很 明显 ， 这 是 一 个 先 增 后 减 的 函数 ， 函 数值 
在 正 半 轴 上 恒 大 于 0， 且 人 数 为 0 时 的 流量 为 0。 

对 拥挤 状态 下 人 员 移 动 速度 与 人 流 密度 的 研究 主要 是 运用 现场 观测 和 录像 
记录 两 种 手段 进行 ,到 目前 已 经 积累 了 大 量 的 观测 数据 。 其 中 ,具有 代表 性 的 研 
究 者 有 前 苏联 的 Predtechenskii Milinskii， 美 国 的 Fruin, Maclennan, Nelson, 
英国 的 Smith, HÆJ Ando， 加 拿 大 的 Paul 等 人 。 各 位 学 者 得 到 的 人 流 密度 与 
行进 速度 关系 存在 较 大 差异 ,他 们 将 人 流 密度 与 行进 速度 的 关系 拟 合成 线性 关系 、 
对 数 关 系 、 指 数 关 系 、 三 角 函 数 关 系 等 。D] 

由 于 要 用 变 分 原理 求 泛 函 极 值 , 且 在 打靶 法 求 边 值 问 题 过 程 中 函数 猜测 值 可 
能 会 超过 函数 在 实际 应 用 中 的 定义 域 , 在 初步 建 模 中 为 了 回避 计算 上 不 必要 的 麻 
Ki, Wi FER BEE SIR LÆR ATA RR OD AEE Ai BERI OS AR e 
同时 ， 为 了 防止 计算 时 提出 的 最 优 解 部 分 超出 实际 应 用 中 的 定义 域 ， 应 有 
p € (一 00,0),v(p) < 0， 如 入流 密度 与 行进 速度 成 负 指 数 关 系 。 但 根据 在 本 校 出 
操 时 对 楼 梯 口 通行 状况 的 观察 ， 这 个 函数 在 p 较 大 时 存在 一 定 失真 。 在 人 员 较 多 
时 ， 楼 梯 仍 然 有 一 个 较 小 的 流通 速度 ， 而 不 是 像 预 测 的 那样 快速 趋向 于 0。 在 本 
校 楼 梯 上 ， 这 个 速度 大 约 为 0.61 人 /s。 

由 此 ， 我 们 想到 演 试 用 截 尾 正 态 分布 的 方法 拟 合 : 设 s(x) 为 待 求 函数 ， 


AT Kar 
s(x) = vo + ve (Sn) — (vo + v,e-CG»*)e7 5 


其 中 ，S 是 通道 面积 ，xi1 为 人 数 ， 其 余 五 个 参数 要 通过 试验 结果 确定 。 第 二 
项 是 罚 函 数 ， 来 保证 人 数 为 0 时 的 流量 为 0。 而 且 ， 第 二 项 的 存在 可 以 使 
p € (795,0), vip) < 0， 保 证 最 优化 解 不 会 超过 函数 在 实际 应 用 中 的 定义 域 。 


3.1.2 实际 测试 


在 对 本 校 楼 梯 的 实际 测算 中 ， 得 到 楼 梯 上 人 数 与 流量 的 关系 如 下 : 


由 于 无 法 发 动 全 校 师 生 为 我 们 做 这 个 试验 , 在 人 数 较 多 时 的 取样 较 少 , 而 人 数 较 
少时 则 获得 了 较 多 样本 。 

拟 合 的 结果 为 : 

k, 20.6519 

k,-5.1020 

p1=1.760 

s=26.77〈 由 测量 得 到 ) 

v0=0.2752 

v1=2.2079 


3.1.3 与 学 界 其 他 结果 的 比较 


Fruin 统计 得 出 下 楼 梯 的 最 佳 条 件 为 :密度 为 2.0 人 /m2, 楼 梯 有 效 宽度 的 流量 
为 1.18 人 /s RI 

中 国 科学 技术 大 学 火灾 科学 国家 重点 实验 室 采用 数码 摄像 的 方法 计算 得 的 
不 同 条 件 下 的 人 流 密度 、 人 流量 、 人 流速 度 参 数值 显示 密度 为 1.8 人 /m2 时 得 楼 
梯 的 最 大 流量 。[3] 

在 人 数 较 多 时 第 二 项 的 影响 可 以 忽略 ， 因 此 : 


Vmax S Vg + Vi = 5.537 


Plv=vmaxr © Pa = 1.760 


与 以 上 两 结果 较为 接近 。 

若 在 人 数 不 过 多 时 以 负 指数 函数 形式 拟 合 速度 与 人 流 密度 关系 ， 同 时 以 
v = plu = ple” -22 形式 估计 人 流量 ， 则 达到 极限 人 流量 时 ， 每 个 人 平均 速度 与 
单个 人 极限 速度 比 满 足 1:e 的 定 值 ， 与 Cu Cs; 无关。 在 本 拟 合 中 ， 这 一 结果 约 为 
1:2.5。 


3.2 条 件 弱 化 模型 


3.2.1 建立 模型 


规定 对 任 一 建筑 物 均 存在 一 特定 时 间 值 t1， 称 为 安全 时 间 ， 在 紧急 情况 发 生 
后 的 时间 内 可 以 保证 建筑 物 的 安全 性 ， 此 后 安全 性 不 能 保证 ， 因 此 我 们 的 目的 
就 是 在 妇 时 间 内 实现 尽 可 能 多 的 疏 散 。 

先 假设 我 们 只 能 控制 每 一 层 的 人 数 ， 不 能 控制 单位 时 间 人 进入 楼 道 的 数量 。 
例如 单位 时 间 进 入 通道 的 人 数 与 时 间 成 负 指 数 函 数 关 系 。 对 于 每 一 段 楼 梯 ， 都 可 
以 对 应 一 个 方程 描述 其 中 人 数 的 变化 。 一 般 的 , 设 人 数 为 x 时 的 疏散 速率 为 s(x)， 
单位 时 间 从 该 层 进入 楼 梯 的 人 数 随时 间 变 化 的 函数 为 y(t)， 则 这 个 方程 的 形式 


为 : 


dx; 
E = —s(x;) + s(xi-4) + yi(t) 
nenne FA A BRUNAR s) = vy + ve PG), þik, 
以 最 简单 的 三 层 为 例 , 将 一 楼 看 作 无 需 经 过 楼 梯 可 直接 逃离 ， 就 可 以 得 到 一 个 方 
程 组 : 


dx, -— 
SEC fi = —S(%4) + ae 
dx» "T 
VU fo = —s(x2) + sa) + (A— a)e 
dx3 
ape 


其 中 表示 已 经 从 通道 出 来 的 人 数 。4 由 总 人 数 决定 ， 而 a 由 分 配方 案 决定 ， 
因此 我 们 可 以 通过 控制 来 找到 最 优 方案 。 


3.2.2 数学 方法 的 证 明 


对 于 一 个 初 值 回 定 的 "元 一 阶 常 微分 方程 组 ， 含 一 个 确定 的 参数 w， 
A: = fi(Xy .Xn T1724 
x; (0) = Xilt=o 
其 图 像 为 n + 12E2E TES [RI ÆR AJI An + 2 维 w， 则 对 于 一 个 确 

定 的 a 将 有 一 条 这 样 的 曲线 ， 所 有 这 些 曲 线 构成 一 个 n + 2 维 空间 中 的 曲面 。 如果 
将 x 也 视 作 变量 是 将 之 前 的 所 有 对 t 求 导 改 为 对 t 求 偏 导 ， 则 形成 边界 条 件 为 在 
t = 0 时 所 有 初 值 已 知 的 偏 微 分 方程 : 

Ox;(t, a) R 

— St) 12 


an 


X;(0,Q) = xi|tzo 
对 于 确定 的 方程 组 ， 其 解 是 唯一 的 。 
形式 地 对 其 求 对 a 的 一 阶 偏 导 ， 得 到 n 个 方程 : 


Oxi(t, a) ` ( Kosi 
EFSER = fi(X1 ...Xn,t,a) i=1,2...n 


由 于 c 不 影响 初 值 ， 所 以 5 = 0 
由 此 可 以 构造 个 辅助 方程 
d ax;(ta) ` mee 
di da 三 万 (0 -Xn t,æ) i = 1,2..n 
Ox; (0,a) 
da id 
Jr an n HO rede e dem + 2 维 空间 中 ， 上 曲面 上 点 在 a 方向 的 全 


X Cl MD. 从 函数 的 角度 看 对 于 每 一 个 组 (t, gq)， 都 有 一 个 确定 的 值 ， 因 此 
可 以 将 系统 看 作对 (t, Qa) 的 二 元 函数 , 对 于 确定 的 时 间 t 求 某 一 个 xi 极 值 的 问题 就 


ECH 0 的 问题 ， 这 是 一 个 两 点 边 值 问题 ， 由 此 可 以 求 出 
x& Aa EL x ARE o 
3.2.3 在 本 模型 中 的 应 用 


根据 偏 微 分 基本 性 质 有 : 


ða Öh 

dt da 

OX 
d'ac Op 


ÒXz 


d da 8 Of 
dt da 
Ox; 
v= (0) =0; T5 (to) = = 0 
RES 
这 需要 运用 到 常 微分 方程 的 数值 运算 和 牛顿 选 代 Qi = Quo, 一 Fe dn 
da? 


REX ， 束 需要 再 形式 求 偏 微 分 一 次 ， 以 求 得 计 


起 更 高 庆 的 问题 ,就 会 出 现 更 多 与 分 村 有 关 的 变量 ,此 时 需要 利用 程 组 的 站 可 
比 矩 阵 进行 牛顿 迭代 : 


A; = Aja —] ME) x E 

有 了 这 个 工具 , 这 个 问题 已 经 基本 上 得 到 了 解决 , 借助 计算 机 软件 可 以 得 到 
足够 好 的 s(x)， 计 算 可 以 选用 打靶 法 或 分 割 法 求 边 值 解 。 

在 本 校 去 年 的 地 震 逃 生 演 习 中 的 观察 表明 , 我 校 教学 楼 四 个 楼 梯 〈 排 布 为 和 矩 
形 四 和 角 ) 中 的 任 一 个 均 有 4 = 10;m = 0.09。 当 时 高 三 部 分 班级 由 于 在 进行 高 考 
复习 没有 参加 演习 ， 因 此 二 、 三 楼 实际 参与 人 数 约 为 450 人 ， 与 按 负 指数 拟 合 
积分 的 结果 大 致 吻合 。 

假设 我 们 可 以 自由 调节 二 、 三 楼 人 数 ， 目 标 时 间 选 取 为 = 60s， 则 通过 以 
上 计算 获得 结果 如 下 : 


Daft 
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80 
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ri] 


我 们 可 以 用 以 上 方法 寻找 qa, 并 对 取 极 值 时 xu xz 对 时 间 的 函数 做 出 精确 的 计 
算 。 但 是 要 画 一 张 精 确 的 xs(to, q) 相 图 的 计算 量 很 大 ， 因 为 这 意味 着 要 对 于 差分 
的 每 一 个 x 值 求 常 微分 方程 的 数值 解 并 绘图 。 此 处 仅 以 步 长 为 1 差分 作 图 ， 给 出 
较为 直观 的 印象 。 图 中 a = 0 代表 所 有 人 分 布 在 2 楼 ，a = 10 代 表 所 有 人 分 布 在 
三 楼 。 以 上 计算 代码 见 求解 代码 @ 

对 于 三 层 以 上 的 建筑 , 同样 可 以 采用 这 一 方法 建 模 并 得 到 各 楼 层 人 数 之 间 的 


建议 比值 ， 唯 一 的 不 同 在 于 方程 数目 以 及 未 知 数 个 数 。 

如 果 要 刻画 人 员 分布 可 变 时 的 下 散 状况 ， 在 之 前 模型 中 有 确定 函数 形式 的 
yb 要 在 新 模型 中 改 为 只 有 确定 边界 值 的 未 知 函 数 。 如 果 我 们 求 出 了 一 个 最 优 玻 
散 方 案 所 要 求 的 y(t) 来 表示 单位 时 间 内 进入 楼 梯 人 数 与 时 间 的 函数 关系 , 我 们 就 
可 以 在 楼 层 中 标定 人 员 跑 到 楼 梯 口 所 需 时 间 t € [ti,t2] 的 对 应 区 间 ， 这 一 区 间 在 
无 险情 时 的 人 员 密 度 应 尽 可 能 接近 上 y(t)dt. 


通过 以 上 数学 模型 与 数学 方法 , 我 们 已 经 可 以 得 到 对 某 建 筑 人 员 分 布 有 一 定 
指导 意义 的 结论 。 在 实际 中 , 我 们 可 以 在 建筑 设计 与 空间 分 配 时 有 意识 的 预 控 每 
一 处 的 人 员 密 度 ， 应 用 此 结论 。 

在 学 校 中 , 由 于 每 个 行政 班 的 人 数 大 致 相同 ， 这 一 模型 的 指导 意义 不 如 之 前 
的 初步 模型 ， 但 是 在 其 他 公共 场合 做 这 样 的 预 控 还 是 可 能 的 。 以 商场 为 例 ,各 店 
面 附近 人 员 密 度 相 差 很 大 , 在 估计 了 每 一 个 楼 梯 所 需要 管辖 的 区 间 与 这 一 区 域 的 
大 致 总 人 数 后 ， 可 以 鼓励 客流 量 不 同 的 商家 选择 与 消防 要 求 相 匹配 的 店面 。 

以 下 尝试 在 y(t) 可 变 的 较 强 条 件 下 ， 求 出 最 优 的 y(t)。 由 于 会 涉及 变 分 法 的 
边界 条 件 ， 以 下 y( 人 9 代表 某 一 时 刻 已 经 进入 楼 梯 的 人 数 ， 单 位 时 间 内 进入 楼 梯 人 
数 与 时 间 的 函数 关系 用 y'(t) 表 示 。 


3.3 条 件 强化 模型 


现在 进一步 考虑 在 控制 每 一 层 人 数 的 同时 , 我 们 也 能 控制 单位 时 间 内 人 进入 
楼 道 的 数量 。 那 么 这 个 问题 就 变 得 复杂 了 ， 需 要 分 多 步 进行 处 理 。 对 于 一 般 的 问 
题 ， 以 三 层 建筑 为 例 ， 玻 散 过 程 可 以 粗略 地 用 这 样 的 方程 组 来 描述 ; 


dr ` 
Dx —s(x4) + 4(t) 


dx; 

kom —s(x2) + (x1) + y«(t) 
ax3 
Pr = s(x2) 


此 模型 可 以 给 我 们 一 个 相对 直观 的 感觉 ， HTH AO yo)» Et 
刻 都 有 一 个 确定 的 xs(to) = 记 " s(x2)dt 与 之 对 应 ,我 们 可 以 定义 一 个 2 = s(x), 


at 
初 值 为 0, 来 代 蔡 积分 形式 。 但 无 论 哪 种 形式 ,都 建立 了 一 个 [Dr(b,》2(b9] 一 xs 的 
了 映射， 由 泛 函 定义 可 知 这 是 二 元 泛 函 ， 我 们 要 求解 的 正 是 泛 函 极 值 。 
这 一 类 中 问题 中 较为 简单 的 问题 如 : 
^E XL BUM y (t) 9 x(to) 


dx li 
as fe y, yt) 
x(0) = xo; YO) = yo; y(t1) = yı 
求 y(t) 使 x(to) 在 其 附近 取 极 值 。 
为 了 解决 这 个 问题 ， 必 须 进行 以 下 几 步 考虑 : 


3.3.1 某 一 瞬间 变化 对 常 微分 方程 的 影响 


EE 常 微分 方程 的 影响 。 定 义 函 数 Fxo,to) 


Z L ft 
x(0) = xo 
f (Xo, to) = x(to) 
o: 
lim f (xo + Aa, e f (Xo, to) 


Aa>0 


为 此 先 讨论 对 于 和 x 等 on ak shea ri El Zoe, 
令 x = x' 十 aq， 原 方程 可 写 为 


+ a) 


=f(x'+a,t) 


其 中 x 初 值 恒 为 xo。 与 之 前 讨论 相同 。 REG: 
ki = d ox ta) dòx' 


dt ða — dt da 
af(x'+a,t) O(x' a) Me ta, t. dr 
有 式 = TE. 
Ox +a da O(x' +a) da 


因此 有 : 


dóx' Of(x rat) Ox 
dt da. a(x’ +a) a 


结合 f xg, to) 定 义 有 : 


m ZO + Aa, e f (xo, to) ELM 
da 


昌 此 一 wem a 与 x 的 影响 可 表示 为 : 
Ax'(ty) = (to) x Aa 


Lu 


Ax(to) = A(x' (t9) + a) = Gm (ty) + 1)Aa 


Ox 
定义 6 = 1, A: 
dp Of(x+Aa,t) 9f(x,t) 


dt  O(x'+Aa) TO B 
B(to) = 
不 妨 把 8 称 为 放大 因子 ， 则 所 时 刻 x 值 的 变化 量 
Ax = BAa 


现在 , 我 们 将 初始 时 刻 微 扰 对 常 微分 方程 的 影响 推广 到 任意 时 刻 微 扰 对 常 微分 方 
程 的 影响 。 定 义 函 数 F(xo, ty, te, 
fut) 
u(0) 2 xo 
x, = u(ti) 


E EG 
dv 
SECH 
v(t,) = uy 
F (Xo, tite, a) = V(L) 
o: 
_ F(X, ty, tè à) — F (xo, ty, t2,0) 
lim ii 一 一 
a>0 a 
fæ: 根据 上 一 步 结 果 ， 设 
QX 
STEIER 
dB  Of(x,t) 


AT A BB) =1 
a > 0 时 对 于 
du' l 
— = fan) 
u'(0) 2 x9 + en 
b(t) 
A: 
u(t) =x ta 
由 唯一 性 定理 
u=v 
THA 
' B(t2) 
v(t;) — u(t;) = u' (t2) — ultz) BŒ) a 


对 于 空 = f (x, 日 在 初 值 zto) = xo + Ad， 其 中 Ac 为 一 阶 小 量 时 ， 这 一 全 
差 是 线性 的 ， 所 以 对 于 tx 时 刻 的 一 个 大 小 为 Aw 的 微 扰 ， 在 所 时 刻 被 放大 为 
ECR Aw， 由 此 可 以 由 一 个 放大 函数 BCb 得 出 任意 时 刻 微 扰 对 最 终 值 的 影响 ， 
在 方程 组 中 计算 这 一 结果 : 


dx ` ; 
a = FYE) 


dy 
Ab go yt) 
仿照 前 面 的 讨论 ， 设 在 如 时 刻 的 两 个 方向 的 微 扰 是 Aaw,Aay， 有 : 
dp, Of af 


dt Re ER 十 ER? 
dYxy g 0g 

dt ES ax Px F ay 7 
dyyx Of 0f 
dt ER + dy PY 


dB, af of 
Cu a Rt m 
其 中 初 值 8 = 1,y = 0。 两 方向 上 的 扰动 : 
Ax = ByAdx + VyxAady 
Ay = es At, + ByAay 
同样 ， 也 可 以 反 过 来 对 确定 


AT 


的 Ax 或 Ay 计 算 对 应 的 Aqx，Aay。 这 说 明 对 于 任 
意 时 刻 的 一 个 微 扰 ， 都 可 以 和 一 个 确定 的 初 值 Ax，Ay 等 交 
代数 的 语言 ， 就 是 : 


Yia 7" Yn 
r -[; . ") 
Yin ^" Ynn 
其 中 yij 代 表 i 方 同 的 微 扰 在 j 方 同 的 影响 。 在 初 值 处 的 矢量 微 扰 被 放大 到 |: 


Yii cU YT Aa, 
i E : |x : 
Vin. Ut" Man Aa, 


对 于 任意 时 刻 的 矢量 微 扰 ， 在 to 时 刻 被 放大 到 : 


Aa, 
T(to) XTX | | 
Aa, 


BOP MARA, Aa,» Ae E AP Fe] EST TE) AAR. DAAE 
在 to 时 刻 的 微 扰 ， 则 之 后 的 方程 就 变 为 了 : 


。 令 Ji = Bi:， 用 线性 


, RAGE 


La — Am) = f(x — Aa, t) 
对 这 个 方程 做 与 之 前 相同 的 等 效 处 理 ， 发 现在 Aai > 0 时 其 仍 会 退化 为 
dB of(x,t) 
dt 0x B 
SCHON, Aa, Aa Z lA AH ASE Bee, RATE TAT Ae TES Ho 
也 可 以 发 现 ， 每 个 微 扰 的 放大 只 取决 于 常量 部 分 ,而 与 小 量 无 关 。 如 果 之 前 的 所 
有 微 扰 所 造成 的 影响 还 是 一 个 一 阶 小 量 的 话 , 它 对 下 一 个 微 扰 所 造成 的 影响 也 没 
有 影响 。 因 此 ， 对 于 一 些 连 续 的 二 阶 小 微 扰 ， 我 们 可 以 对 它们 进行 积分 。 


3.3.2 连续 二 阶 小 量 对 常 微分 方程 的 影响 


可 以 注意 到 , 这 里 的 二 阶 小 微 扰 由 y 的 变 分 提供 , 积分 后 得 到 的 是 一 阶 小 量 
不 影响 之 后 的 扰动 。 不 妨 先 假设 这 个 二 阶 小 的 微 扰 是 显 式 ， 很 容易 得 到 : 


Ey 


MA 
Ax = Des A? a(t)dt 


d d GH GE a 
FII = f(x,y, Yy , D EIN Ay (MER eo, 其 在 dt 对 应 的 瞬间 微 扰 为 


af af 


Cy by + 5,:8y dt 
由 于 系统 中 显 含 x 这 一 微 扰 最 终 被 放大 到 | 
BO) òf Of 


BO e 十 po ER 


up) of Of |, 
|, pq ay? * Gyro Xt 


np) of Of |, 
, eo ay” * ay 7 dt 


5 B(t,) af TIO 
-Ja BO ay * J, pO ay” A 


pup) of BD OF D, d (BD of 
=J. Bo ay trO ;"» Å oy x Za 


(t1) 0f 
= 6y(t1) = 6y(t2,* EET yl = 0 


BW) dy dt PE) dy’ 


— Sy + 
. Bc) ay?” * ay! 
~ Oyit 


2 [2D ai Baue [^ Gou d BC) 2) SS 


BO dy dt BO ay’ ` 
HEAD (t) FT DAT, 但 由 于 在 多 变量 时 转 置 矩阵 是 不 可 约 的 ， 所 以 先 写成 
这 样 保证 形式 。 整 理 一 下 ， 对 于 一 个 确定 的 y(t》， 我 们 有 : 


SE EE n 


dx ` A 
a Í, 
dB OfCQoyyt) B 


IK Tar bi SET 


, B(t)dy dt (t) ay’ 
当 它 取 极 值 时 ， 有 : 


BA d f, of _ 
Bo dtp oy 


从 中 我 们 可 以 获得 2 的 显 式 ， 约 东 条 件 为 ; 


X(to) = xo, y(to) = Yo, Y(&) = yo P (to) = 1, Bi = P (ti) 
这 样 就 再 次 形成 常 微分 方程 的 两 点 边 值 问题 。 


3.3.3 建 模 过 程 


由 3.1 讨论 知 : 
"ean —k4 (p1)? EE 
mente — (vo + vye ON )e 5 
VER BE Er ESTE T8] PI AGENT EE AGE, EE NFT) ESÈ JA 
li vox,e dt DAR de ME ERT, RARA i OT AKI te BBE 


内 。 否 则 我 们 大 可 以 在 t = 0 时 刻 让 一 大 批 人 在 尽 可 能 短 的 时 间 内 进入 楼 梯 ， 以 
使 楼 梯 的 朴 散 速率 达到 极 大 值 ， 而 完全 不 顾 这 一 举动 对 楼 梯 口 所 造成 的 危险 ,这 
就 是 3.4 开始 的 模型 符合 数学 工具 的 要 求 但 实际 上 并 不 合理 的 原因 。 

大 多 数 公共 场合 的 楼 梯 都 设 有 防火 门 、 转 角 等 障碍 ,即使 是 与 楼 层 直 接 相连 
的 楼 梯 , 在 进入 时 人 流 宽度 也 会 减 小 一 半 , 因此 此 处 也 成 为 了 各 种 意外 的 高 发 地 。 
在 定量 估计 其 造成 的 风险 时 ， 可 以 视 风 险 Rer 为 

Rer = A(Ns/Ne) Fy 
《 取 自 《人 和 群 拥挤 踩踏 事故 后 果 微 观 建 模 及 模拟 分 析 》 式 中 Ns 为 + 时 刻 特定 区 域 
内 的 滞留 人 数 , Fy 为 t 时 刻 特定 区 域 拥挤 人 群 中 的 死亡 人 数 , Ni 为 设计 人 数 , AN 
触发 因子 ， 可 通过 专家 打分 和 评定 得 到 ) 5 

在 本 模型 中 , 辕 于 计算 能 力 和 调查 能 力 , 近似 设 单位 时 间 的 涌 入 人 数 正 比 于 
这 段 时 间 内 入 口 附近 人 数 ， 即 


y«N, 

原文 中 的 是 由 拥挤 力 达 到 242N 持续 3min 时 造成 人 员 伤 亡 判 定 的 ， 而 本 模 
型 中 假设 伤亡 以 外 的 因素 也 会 导致 判定 ， 因 此 取消 3min 的 延迟 时 间 〈 这 一 假设 
也 会 带 来 误差 , 因为 长 期 拥挤 带 来 的 体力 下 降 也 可 能 对 意外 发 生 的 期 望 有 影响 )， 
只 认为 X Ns。 

因此 , 任意 时 刻 的 风险 预期 可 表示 为 41(y)2, 其 中 参数 待定 。 对 于 本 校 楼 梯 ， 
可 以 认为 4 = 0.04〈 这 一 结果 是 通过 观察 本 校 学 生 第 三 节 课 下 课时 下 楼 抢购 课 间 
餐 的 行为 进行 估计 的 ,不 精确 之 处 难以 避免 ,希望 在 有 条 件 时 能 得 到 更 精确 的 结 
果 使 得 模型 的 实际 意义 更 强 )。 这 个 量 对 时 间 积 分 就 是 全 下 散 过 程 中 总 风险 的 期 
望 。 

对 一 个 玻 散 方案 作 评 价 时 应 同时 考虑 到 通道 的 情况 和 入 口 的 风险 , Box AI VE 


£ Ll) E 
价 函 数 ， Mix, = fe Vo + Vie (s pi) — (Yo + ve ™ PD le 


5-20). 5 
成 微分 形式 ， 则 ; 


dx sl) -kaxi 
E = —vy — vie (Son) 二 (vo + ve ED es ty 
dx» 


2 k 
= Vo 十 vye fè-n.) — (vo + vye ka ye — A'Y” 
x,(0) = x,(0) = 0 
y(0) = 0, y(to) = Yo 


dr 


求 y(t) 使 x， (to) 取 极 值 。 
以 下 以 s= 10;v0 = 0.182;v1 = 1.46; 的 小 通道 , 管辖 半径 内 最 远 的 人 在 50s 内 到 达 
AH, f 80 人 为 例 计 算 。 


3.3.4 数学 处 理 以 及 在 matlab 平台 的 实现 


a] 3.4 的 讨论 ， 对 于 一 个 初始 时 刻 的 微 扰 


Fe 
Aa» 
会 被 放大 到 
E 
因此 在 t 时 刻 
Fa Gi 
Aa, (t) 


可 以 等 效 为 t= 0 时 刻 的 
人 T ( m x a 


Aa, (0) Ya2 V22 Aa, (t) 
其 变 分 为 : 

d d 

ti -1 Of sy 4 A syr 

| Ge e AE a » du du dt 

o Wiz(ti) Y22 (ty) Yi2 V22 Oo PE 
dy y ðy' y 
dëi 


E e a 2 du 
Yaz(03) Y22 (ti) Cer, V22 


5 
i=1 == 


ii 
0 


d | (vii(t1i) A Vi1 Ya du 
me Ou; dt 
dt br V22(t1) 2 Es EJ a Ui 


由 于 6ui 任 意 性 ， 为 使 对 于 任意 x2(ti) 取 极 值 ， 有 : 


Gow e SH dy 


x 
Ya2() Y22(ti) Vi2 V22 9f; 
Qy 
an 
a Var. ee SN du 
dt| WizCti) Y22(ti) Vin. V22 af, 
dy’ 


的 第 二 行为 0 
计算 出 满足 边界 条 件 的 解 如 图 所 示 


在 计算 的 开始 部 分 图 像 有 大 范围 震荡 , 这 是 由 于 计算 方法 引起 误差 。 在 修改 


初始 猜测 值 时 震荡 的 大 小 和 方向 都 随 之 变化 ,而 后 面 主 要 的 图 像 不 变 。 其 中 第 三 
张 图 像 对 某 一 通道 所 管辖 半径 内 人 员 的 分 布 有 一 定 指导 意义 。 

和 其 他 以 迭代 为 基础 的 方法 一 样 , 本 方法 的 成 功 与 否 很 大 程度 上 依赖 于 初 值 
是 否 能 使 迭代 向 所 需 方向 收敛 , 而 且 由 于 需要 猜测 的 包括 转 置 矩阵 的 元 素 , 这 些 
元 素 在 y( 提 未知 时 较 难 估计 ， 因 此 对 初 值 猜测 的 要 求 较 一 般 欠 代 问题 要 高 。 

如 果 初 值 选 取 不 当 ， 也 会 出 现 这 种 解 : 


dydt 


dy/dt 


以 上 计算 代码 见 求 解 代码 @ 


本 校 实际 玻 散 情 况 为 有 4 个 对 称 分 布 在 矩形 四 角 的 4 个 楼 梯 。 易 证 优化 方 
案 一 定 满足 四 个 楼 梯 之 间 的 对 称 。 假 设 2、3 楼 之 间 的 楼 梯 压 力 较 小 ， 暂 时 优先 
考虑 1、2 楼 之 间 楼 梯 承 受 的 压力 ， 这 样 ， 我 们 需要 求 的 2 楼 之 间 楼 梯 单 位 
时 间 进 入 人 数 与 时 间 的 关系 同样 可 以 由 以 上 方程 组 定义 的 泛 函 极 值得 到 。 代 入 之 
前 测 得 的 数据 : 
k1=0.6519 
k,=5.1020 
p1=1.760 
s=26.77( 由 测量 得 到 ) 
v0=0.2752 
v1=2.2079 
解 的 图 像 如 下 : 


bu 


140 


- 100 


10) 


dydt 


dyidt 


以 上 计算 代码 见 求解 代码 @ 

如 果 认 为 2、3 层 间 的 楼 梯 所 承受 压力 也 要 计 入 ， 则 需要 再 添加 方程 。 事 实 
上 ， 和 3.2 的 模型 一 样 ， 对 于 3 层 或 以 上 的 朴 散 同样 可 以 利用 变 分 原理 得 到 解 ， 
只 需要 添加 方程 与 未 知 数 即 可 , 但 是 在 解 的 过 程 中 要 引入 泛 函 极 值 的 变 边 值 方法 。 
如 三 层 问 题 的 所 有 方程 与 边界 条 件 为 : 


k2x1 


dx ET. S - 
2 = —Vo _ vie KU pi) + (Yo + ve k(x)" le S + yu 
dx -k (22-0) -kaxa pp ) 
T = -v — we ki (5 pi) 十 (vo a We iby le S +Vo+14€ ki (5 pi) 
|k2xX4 
= (Yo + ve 100" Je S + Va 
dx -k (X20 JEG 
yt me a($-A) — (v, + ve MP) — A7? + 02) 


x1(0) = x2(0) = x3(0) = 0 
yi (0) = ya (0) = 0; y1 (to) + Y2(to) = yo 
Ky. yo xs (to) BUE. ES 3 x 3 的 转 置 矩阵 [的 方程 ， 方 法 如 3.4.1, 
这 里 不 再 展开 。) 对 应 的 欧 拉 方 程 组 


JA of 
Oy; Oy; 
òf, d òf, 
= Tes mM —1 NL) i 二 
T(to) x TCE x Ze 3 T(to) x IT (t) x ay, ,i= 1,2 
dë dë 


Oy; Oy; 


的 第 三 行为 0， 因 此 又 增加 了 两 个 方程 ， mam = = yi, i = 1,2 以 使 其 成 为 动 


力 系 统 ， n. 16 个 方程 ， 


有 原始 边界 条 件 6 个 , 转 置 矩 阵 初 值 等 于 单位 矩阵 与 转 置 矩阵 终 值 等 于 猜测 矩阵 


同时 转 置 矩阵 的 终 值 共 9 个 元 素 为 待定 参数 。 方 程 中 


各 页 献 9 man 3 共 24 个 ， 还 缺少 一 个 。 这 一 个 是 由 泛 函 极 值 的 变 边 值 方 


MLB]. MARI 


OF (to) = T (to) x I! (t) 


法 决定 的 ， 不 难 发 现 由 于 限制 条 件 yj (to) + yz(to) = yo， 两 个 端点 变 分 并 不 是 独 


这 个 矩阵 的 第 三 行为 0。 


$ 


T(to) Xf (f) 


结 Ay (to) + yz(to) = = yo» 得 到 最 后 一 个 边界 条 件 : 


0f, 0f, 
Oy dy,’ 
Òf, —1 af, 
ER y1 Cto) + T(to) x (t) ovs by2(to) 
dr Of 
dy,’ dy,’ 

of 0f, 

Oy Oy2 

af, -1 fo 

TA Au (to) — T (to) x (t) By; 6y2(to) 

dr Of 

Oy dy," 


这 个 矩阵 的 第 三 行为 0。 


由 此 可 见 ,通过 变 分 | 


在 变量 和 函数 增多 时 差分 法 所 需 差分 的 维度 增高 ， 误 差 将 呈 几 何 级 数 放 大 ,这样 


原理 求 得 优化 方案 在 理论 上 是 可 行 的 .和 3.2 模 型 一 样 ， 


的 方法 可 以 弥补 这 一 缺陷 。 虽 然 计 算 量 较 大 , 但 是 它 比 直接 差分 法 寻找 最 优点 ( 寻 
找 函 数 时 整个 函数 的 定义 域 都 需要 差分 ) 或 无 目的 尝试 要 简便 许多 。 


总 结 


本 文采 用 分 析 的 方法 研究 这 一 模型 , 在 大 部 分 情况 下 收 到 了 较 好 的 效果 , 但 
是 和 所 有 采用 迭代 方法 的 计算 一 样 ， 这 一 方法 也 依赖 于 选取 初 值 的 收敛 性 , 所幸 
的 是 一 般 的 足 散 模型 中 函数 的 凹凸 性 都 较 单 一 ， 否 则 和 帮 代 法 的 成 功率 会 较 小 。 

本 文 的 模型 建立 关于 时 间 和 研究 条 件 的 限制 , 采集 的 数据 量 不 够 大 , 但 是 对 
楼 梯 琉 散 问 题 的 研究 该 方法 具有 普 适 性 , 如 果 能 将 其 与 足够 精确 的 预报 模型 结合 ， 
或 将 使 得 模型 具有 更 大 的 实践 价值 。 
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附录 


求解 代码 : 
(03 层 建筑 中 2、3 MEN Å DIA LG OA BLE 


function r = bf3 


kl = 0.6519; 
roul = 1.7600; 
s — 26.7728; 
v0 = 0.2752; 
vl = 2.2079; 


tO = 60; 
m = 0.09; 
A = 10; 


lambda = 0:.1:1; 
a = A*lambda; 


function x3t0 = ff (alpha) 


sfun = @(x) vO + vl.*exp(-kl.*(x./S - roul).^2); 


figure (1); 
plot(0:100,sfun(0:100),'k:'); 
xlabel('S$\mathrm{x}$','fontsize',15, 'interpreter','latex'); 


ylabel('S\mathrm{s}$','fontsize',15, 'interpreter','latex'); 


title ('S\mathrm{s}{(x)}$','interpreter', 'latex', 'fontsize',15); 


function dx1 = DxiDtSubFun(t, x1) 
dx1 = -sfun (x1) + alpha.*exp(-m*t); 


end 


tspan = 0:60; 

%tspan = [0 60]; 

options = odeset('RelTol', 0.0001); 

ode45 (@Dx1DtSubFun, tspan, 0, options); 


sol 


t45 


sol.x; 


x45 = sol.y; 


figure(2);clf; hl = axes; %hold on; 
plot (hl, t45, x45, 'r:'); 
xlabel("Nitt'); ylabel('\itx 1'); title('$x 1(t)$', 


'interpreter', 'latex', 'fontsize', 15); 


h = .1;%2%3n 
tneed = 0:h:t0; 


[valxl,derx1] = deval(sol, tneed); 


= derxl + A*exp(-m*tneed) ; 


Q 


x2 = zeros(length(tneed),1); 
for i = 1:length(tneed)-1 


Kl = h*(-sfun(x2(i)) + g(i)); 

K2 = h* (-sfun(x2(i)+ K1/2) + (g(i)+g(i+1))/2 ); 

K3 = h* (-sfun(x2(i)+ K2/2) + (g(i)+g(i+1))/2 ); 

K4 = h*(-sfun(x2(i)+ K3) + g(i+l)); 

x2(i+1) = x2(i) + 1/6*(K1 + 2*K2 + 2*K3 + K4); 
end 


figure (3); 
plot(tneed, x2, 'r:'); 
xlabel("Nitt'); ylabel('\itx 2'); title('Sx 2(t)S', 


'interpreter', 'latex', 'fontsize', 15); 


end 


$»y'OCÓx3|tO0 

figure(5) 

plot(tneed, sfun(x2)); 

x3t0 = h*sum( sfun(x2(1l:end-1))+0.5*diff (sfun (x2) ) 


x3 = arrayfun(@ff,a); 


figure (4); 


plot(a, 


xlabel ( 


x3, tk) 


'Nita'); ylabel('Nitx 3'); title('$x 3(t0,a)S', 


'interpreter','latex', 'fontsize', 15); 


end 


GmbH ITZ PR TEE BE) ORG rh DA A 


function r = bf 


lambda 


= 0.015; 


)4 


kl = 0.6519; 
k = 5.1020; 
roul - 1.76; 
s = 26.7728; 
0.61375 
4.9234; 


vO 


< 
pa 
ll 


tO = 70; 
vt0 = 300; 
step = 0 2 


x10 = 0; 

x20): = OY 

yO = 0; 

ml = 2*kl*vl/s; 

m2 = k/s*(vO*vl*exp(-kl*roul^2)); 


$Df = dfl/dx1l£' E«p£O 

Df = Q(xl)ml.*(x1l./s-roul).*exp(-kl.*(x1l./s-roul1).^2) - 
m2.*exp(-k./s.*x1); 

figure (3); 

plot(0:100,Df(0:100),'k:');xlabel('x1'); ylabel('Df'); 
title('df1/dx1'); 


oO 
oO 


function dx1 = Dx1DtSubFun(t, x1) 
dx1 = 
-v0-vl*exp(-kl*(xl/s-roul).”2)+(vO0+vl*exp(-kl*roul”2))*exp(-k*x1/s); 
dy = zeros(size(t)); 
if length(t)>1 
for j = I:length(t) 
if SIT € I 


dy(j) = -(beta(t0O)*exp(-t(j)*Df(0.5*t(3))) + c2) / 
(2* lambda) ; 
else 
dy(j) = -(beta(t0) /beta(t(j))+c2)/(2*lambda); 


end 


end 


TE t « 1 
dy = -(beta(t0)*exp(-t*Df(0.5*t)) + c2) / 
else 
$ if isnan(beta(t)) Sé -isnan(beta(t0)) 
% disp(['t:',num2str(t)]); 
% disp(['beta(t):',num2str(beta(t)),' 
beta(t0):',num2str(beta(t0))]); 
$ end 
dy = -(beta(t0) /beta(t)tc2) /(2*lambda); 
end 
end 


dx1 = dx1 + dy; 


end 


epsilon = 10^(-3); 
beta = @(t)exp(-0.003*t) ;%?A2426ER%Op-0.0001 på ú5 I 
betafirst = beta; 


tspan = [0, t0+10]; 
x10 = 0; 


for kk = 1:2 
for k = 1:Num 


tao = O:step:t0; 
tmp = 0; 
for i = I:length(tao)-l 


if tao(i)«121*ÁEEY?ywarning 


tmp tmp + step*exp(-tao(i)*Df(0.5*tao(i))); 


else 


tmp tmp + step*1/beta(tao(i)); 


end 
end 


integ = tmp;%>ý Ö 


c2 = (-beta(t0) * integ - 2*lambda*yt0 )/t0; 
$c2 = o2 9 


dy0 = -1/ (2*1ambda)* (beta (t0) +c2); 
dyt0 = -(1tc2) / (2*1ambda); 
dx1t0 = dyt0 - lambda*dyt0^2; 


(2* lambda) ; 


options = odeset('RelTol', 0.001); 

sol = ode45(@Dx1DtSubFun, tspan, x10, options); 
t45 = sol.x; 

x45 = sol.y; 


figure(1);clf; hl = axes; Shold on; 

figure(2);clf; gl = axes; Shold on; 

plot (hl, t45, x45, 'r:'); xlabel(hl, 't'); ylabel(hi, 'xl'); title(hl, 
"xl (t)'); 

plot(gl, x45, Dx1DtSubFun(t45, x45), 'b-'); xlabel(gl, 'Í»OÓEx1'); 
ylabel (g1, 'FÜqEdx1'); title(gl, '"xluÄlamsAæ"):; 


[x45t0,dx45t0] = deval(sol, t0); 


dy = G(tt)-(beta(tO)/beta(tt) + c2) / (2*lambda); 
z = O:step:t0; 

y = zeros (size(z)); 

dydt = zeros(size(z)); 


for i = 1:length(z) 


if z(i)«.591?ÁEEY?ywaring 
dydt (i) -(beta(t0)*exp(-z(i)*Df(0.5*z(i))) + c2) / 
(2* lambda) ; 
else 
dydt (i) 


dy(z(i)); 


end 


vill = 0; 


continue; 


y(i) = y(i-1)+ step*0.5* (dydt (i) +dydt (i-1));%»y-0 


(z, y, 'b:');xlabel('t'); ylabel('y'); title('y(t)'); 


plot(z, dydt, 'b:');xlabel('t'); ylabel('dy/dt'); title('dydt'); 


dydtsquare 
x2tOtmp = -x45t0 + yt 
+ 0.5*diff(dydtsquare)) 


if k == 
x2t0 = x2tOtmp; 
bestnum = k; 
end 


dydt.^2; 


0 - lambda * step * sum(dydtsquare (1:end-1) 


if x2t0 < x2tOtmp & x2t0tmp > 0 


x2t0 = x2t0tmp; 
bestnum = k; 
end 
end 
betatmp = @(t)Df(deval(sol, t)); 
beta = @(tt)exp(quadl(betatmp, 0, tt)); 
end 
if kk == 
Num = bestnum; 
disp(['bestnum:',num2str (bestnum)]); 
beta = betafirst;$?À?&5*6fÉWwOÓpn 
end 
end 
end 


@ 条 件 较 强 模 型 的 泛 函 方法 在 本 校 教 学 楼 消防 中 的 应 用 尝试 


function r = bf2 
lambda = 0.04; 
kl = 0.6519; 

k = 5.1020; 

roul - 1.76; 

S = 26.7728; 

v0 = 0.2752; 

vl = 2.2079; 


tO = 70; 
vt0 = 191; 
step = 0 9 


x10 = 0; 
x20 = 0; 
yO = 0; 


ml = 2*kl*vl/s; 
k/s*(vO0+vl*exp(-kl*roul”2)); 


可 
N 
ll 


$Df = dfl/dxl£ E«u4£06 

Df = QG(xl)ml.*(x1l./s-roul).*exp(-kl.*(x1l./s-roul1).^2) - 
m2.*exp(-k./s.*x1); 

figure (3); 

plot (0:100,D£(0:100),'k:');xlabel('x1'); ylabel('Df'); 
title ('afl/axl'); 


ob 
ER 


function dx1 = Dx1DtSubFun(t, x1) 
dx1 = 
-v0-vl*exp(-kl*(xl/s-roul).”2)+(vO0+vl*exp(-kl*roul”2))*exp(-k*x1/s); 
dy = zeros(size(t)); 
if length(t)>1 
for j = I:length(t) 
TE EG) < 1 


dy (j) = -(beta(t0O)*exp(-t(j)*Df(0.5*t(3))) + c2) / 
(2* lambda) ; 
else 
dy(j) = -(beta(t0) /beta(t(j))+c2)/(2*lambda); 
end 
end 
else 
TÉ t «I 
dy = -(beta(t0)*exp(-t*Df(0.5*t)) + c2) / (2*lambda); 
else 


ob 


if isnan(beta(t)) && -isnan(beta(t0)) 
disp(['t:',num2str(t)]); 


ob 


% disp(['beta(t):',num2str(beta(t)),' 


beta(t0):',num2str(beta(t0))]); 
$ end 
dy = -(beta(t0) /beta(t)+c2)/(2*lambda); 
end 
end 


dx1 = dx1 + dy; 


end 


epsilon = 10^(-3); 
beta = Q(t)exp(-0.003+t);$32A?A3GÉkMÒp-0.0001 på ú5 I 
betafirst = beta; 


tspan = [0, t0+10]; 
x10 = 0; 


for kk = 1:2 
for k = 1:Num 


tao = O:step:t0; 
tmp = 0; 
for i = 1:length(tao)-1 


if tao(i)«121*ÁEEY?ywarning 


tmp tmp + step*exp(-tao(i)*Df(0.5*tao(i))); 


else 


tmp tmp + step*1/beta(tao(i)); 
end 
end 


integ = tmp;$»y-O 


c2 = (-beta(t0) * integ - 2*lambda*yt0 )/t0; 
$c2 = c2 — 1 


dy0 -1/ (2* lambda) * (beta (t0)+c2); 
dyt0 = -(1tc2) /(2*lambda); 
dx1t0 = dytO - lambda*dyt0^2; 


options = odeset('RelTol', 0.001); 
ode45(8DxlDtSubFun, tspan, x10, options); 


sol 
t45 


sol.x; 


x45 = sol.y; 


FIG 
fig 
plo 
'x1 (t) 
plo 
ylabel 


[x4 


Lo 
ll 


dyd 


for 


ure(1);clf; hl = axes; %hold on; 
ure(2);clf; gl = axes; Shold on; 
t(hl, t45, x45, 'r:'); xlabel(hl, 't'); ylabel(hi, 'xl'); title(hi, 


D. 
t(gl, x45, DxiDtSubFun(t45, x45), 'b-'); xlabel(gl, 'Í»OEx1'); 
(gl, 'ÉÒfÈdxl'); title(gl, 'xlpAIaFAe') ; 


5t0,dx45t0] = deval(sol, t0); 


G(tt)-(beta(tO)/beta(tt) + c2) / (2*lambda); 
0:step:t0; 

zeros(size(z)); 

t = zeros(size(z)); 


i = 1:length(z) 


if z(i)«.591?ÁEEY?ywaring 
dydt (i) -(beta(t0)*exp(-z(i)*Df(0.5*z(i))) + c2) / 


(2* lambda) ; 


else 
dydt (i) 


dy (z (3)); 


end 
if i== 
vill = 0; 


continue; 


y(i) = y(i-1)* step*0.5* (dydt (i) +dydt (i-1));$»y-Ó 


(z, y, 'b:');xlabel('t'); ylabel('y'); title('y(t)'); 


t(z, dydt, 'b:');xlabel('t'); ylabel('dy/dt'); title('dydt'); 


dydtsquare = dydt.^2; 
x2tOtmp = -x45t0 + ytO - lambda * step * sum(dydtsquare (1:end-1) 


+ 0.5*diff(dydtsquare)) 


x2t0 = x2tOtmp; 
bestnum = k; 


end 


if x2t0 < x2tOtmp & x2t0tmp > 0 
x2t0 = x2t0tmp; 
bestnum = k; 


end 


betatmp = @(t)Df(deval(sol, t)); 


beta = @(tt)exp(quadl(betatmp, 0, tt)); 
end 


Num = bestnum; 
disp(['bestnum:',num2str (bestnum)]); 


( 
beta = betafirst;$?À?8?6ÉwÓgu 


end 


end 


end 


